#!/usr/bin/env python 
# -*- coding: utf-8 -*-
# @Project : model
# @Author  : 邢文成
# @File    : solve.py
# @Software: PyCharm

import copy
import matplotlib.pyplot as plt
import numpy as np
import xlrd
import math

from question1.NSGA2 import NSGA2

pa = None
od = None
R = None
L = None
se = None
cw = 1
cd = 1
a = 1
b = 1
result = set()


class Gene(object):

    def __init__(self, m, f, hs, he):
        self.m = m
        self.f = f
        self.theta = 1 / (m + 1)
        self.hs = hs
        self.he = he
        self.CP = None
        self.CO = None
        self.yes = False

    def fit(self):
        if self.yes:
            self.CP = self.pro_CP()
            self.CO = self.pro_CO()
        else:
            self = generate()
            self.CP = self.pro_CP()
            self.CO = self.pro_CO()

    def constraint(self):
        z = 1860
        w = np.zeros(31, dtype=np.int_)
        for k in range(1, 31):
            if self.hs <= k <= self.he:
                w[k] = 0
        tmp = []
        for k in range(1, 31):
            for n in range(1, k + 1):
                for j in range(k + 1, 31):
                    tmp.append(w[k] / z * od[n][j])
        q1 = max(tmp)
        tmp1 = []
        for k in range(1, 31):
            for n in range(self.hs, k + 1):
                # print(self.hs)
                # print(self.he)
                # print(k + 1)
                for j in range(k + 1, self.he + 1):
                    tmp1.append((1 - self.theta) / z * od[n][j])
        # print(tmp1)
        q2 = max(tmp1)
        T = 3600
        Iz = 108
        if self.f >= q1:
            if self.m * self.f >= q2:
                if self.m * self.f + self.f <= T / Iz:
                    tmp_list = []
                    for i in range(1, 31):
                        if 1 <= i <= self.hs or self.he <= i <= 30:
                            # print(pa[i][i+1])
                            tmp_list.append(pa[i][i + 1])
                            Bi = self.pro_Bk(i)
                            tmp_list.append(Bi)
                            Di = self.pro_Dk(i)
                            tmp_list.append(Di)
                            # if max(Bi, Di) / z <= self.f <= max(tmp_list) / z:
                            # print('# if max(Bi, Di) / z <= self.f <= max(tmp_list) / z:')
                    # print(self.f)
                    # print(max(tmp_list))
                    if self.f >= max(tmp_list) / z:
                        self.yes = True
        else:
            self.yes = False

    def __str__(self):
        return "m:{}\nf:{}\nhs:{}\nhe:{}\nCP:{}\nCO:{}\n".format(self.m, self.f, self.hs, self.he, self.CP, self.CO)

    def __hash__(self):
        return hash("".join(map(str, str(self.m) + str(self.f) + str(self.hs) + str(self.he))))

    def __eq__(self, other):
        if self.m == other.m and self.f == other.f and self.hs == other.hs and self.he == other.he:
            return True
        else:
            return False

    def copy(self):
        indiv = Gene(self.m, self.f, self.hs, self.he)
        indiv.CP = copy.deepcopy(self.CP)
        indiv.CO = copy.deepcopy(self.CO)
        indiv.yes = copy.deepcopy(self.yes)
        return indiv

    def pro_CP(self):
        CW = self.pro_CW()
        CVR = self.pro_CVR()
        CVD = self.pro_CVD()
        return cw * CW + cd * (CVR + CVD)

    def pro_CW(self):
        CW1 = self.pro_CW1()
        CW2 = self.pro_CW2()
        return CW1 + CW2

    def pro_CW1(self):
        T = 3600
        PA1 = self.pro_PA1()
        return T / (2 * self.f) * PA1

    def pro_CW2(self):
        T = 3600
        PA2 = self.pro_PA2(hs=self.hs, he=self.he)
        # print(PA2)
        # print(self.m)
        return (T * (self.theta + (1 - self.theta) / self.m)) / (2 * self.f) * PA2

    def pro_PA1(self):
        PA2 = self.pro_PA2(hs=self.hs, he=self.he)
        sum_s = 0
        for k in range(1, 30):
            for j in range(k + 1, 31):
                sum_s += od[k][j]
        return sum_s - PA2

    def pro_PA2(self, hs, he):
        sum_s = 0
        for k in (hs, he):
            for j in range(k + 1, he + 1):
                sum_s += od[k][j]
        return sum_s

    def pro_CVR(self):
        sum_s = 0
        for k in range(1, 31):
            for j in range(2, 31):
                sum_s += od[k][j] * R[k][j]
        return sum_s

    def pro_CVD(self):
        CVD1 = self.pro_CVD1(hs=self.hs)
        CVD2 = self.pro_CVD2(hs=self.hs, he=self.he)
        CVD3 = self.pro_CVD3(hs=self.hs, he=self.he)
        CVD4 = self.pro_CVD4(hs=self.hs, he=self.he)
        CVD5 = self.pro_CVD5(hs=self.hs, he=self.he)
        CVD6 = self.pro_CVD6(he=self.he)
        return CVD1 + CVD2 + CVD3 + CVD4 + CVD5 + CVD6

    def pro_CVD1(self, hs):
        sigma = 0.04
        tmp = 0
        for k in range(1, hs):
            for j in range(k + 1, hs + 1):
                tmp1 = 0
                for n in range(k, j):
                    Dn = self.pro_Dk(n)
                    Bn = self.pro_Bk(n)
                    tmp1 += (Dn + Bn) / self.f
                tmp = od[k][j] * tmp1
        return sigma * tmp

    def pro_CVD2(self, hs, he):
        sigma = 0.04
        tmp = 0
        for k in range(1, hs):
            for j in range(hs + 1, he + 1):
                partL = 0
                partR = 0
                tmp1 = 0
                for n in range(k + 1, hs):
                    Dn = self.pro_Dk(n)
                    Bn = self.pro_Bk(n)
                    tmp1 += (Dn + Bn) / self.f
                partL = tmp1
                tmp2 = 0
                for n in range(hs, j):
                    for q in range(hs + 1, j + 1):
                        tmp2 += 2 * od[n][q] * 1 / (self.f + self.m * self.f)
                partR = tmp2
                tmp += od[k][j] * (partR + partL)
        return sigma * tmp

    def pro_CVD3(self, hs, he):
        sigma = 0.04
        sum_s = 0
        for k in range(1, hs):
            for j in range(he + 1, 31):
                tmp1 = 0
                partL = 0
                partM = 0
                partR = 0
                for n in range(k, hs):
                    Dn = self.pro_Dk(n)
                partL += Dn / self.f
                for n in range(hs, he):
                    for q in range(n + 1, he + 1):
                        partM += 2 * od[n][q] * 1 / (self.f + self.m * self.f)
                for n in range(he, j):
                    Dn = self.pro_Dk(n)
                    partR += Dn / self.f
                sum_s += od[k][j] * (partL + partM + partR)
        return sigma * sum_s

    def pro_CVD4(self, hs, he):
        sigma = 0.04
        sum_s = 0
        for k in range(hs, he):
            for j in range(k + 1, he + 1):
                tmp = 0
                for n in range(k, j):
                    Dn = self.pro_Dk(n)
                    Bn = self.pro_Bk(n)
                    tmp += (Dn + Bn) * 1 / (self.f + self.m * self.f)
                sum_s += tmp * od[k][j]
        return sigma * sum_s

    def pro_CVD5(self, hs, he):
        sigma = 0.04
        sum_s = 0
        for k in range(hs, he):
            for j in range(he + 1, 31):
                tmp = 0
                for n in range(k, j):
                    Dn = self.pro_Dk(n)
                    Bn = self.pro_Bk(n)
                    tmp += (Dn + Bn) / self.f
                sum_s += od[k][j] * tmp
        return sigma * sum_s

    def pro_CVD6(self, he):
        sigma = 0.04
        sum_s = 0
        for k in range(he, 31):
            for j in range(k + 1, 31):
                tmp = 0
                for n in range(k, j):
                    Dn = self.pro_Dk(n)
                    Bn = self.pro_Bk(n)
                    tmp += (Dn + Bn) / self.f
                sum_s += od[k][j] * tmp
        return sigma * sum_s

    def pro_Dk(self, k):
        sum_s = 0
        for j in range(k + 1, 31):
            sum_s += od[k][j]
        return sum_s

    def pro_Bk(self, k):
        sum_s = 0
        for j in range(1, k):
            sum_s += od[j][k]
        return sum_s

    def pro_CO(self):
        sum_s = 0
        for h in range(1, 3):
            Kh = self.pro_Kh(h)
            sum_s += Kh
        return a * (1 + self.m) * self.f + b * sum_s

    def pro_Kh(self, h):
        if h == 1:
            return self.pro_K1()
        elif h == 2:
            return self.pro_K2(hs=self.hs, he=self.he)

    def pro_K1(self):
        sum_s = 0
        for i in range(1, 30):
            sum_s += L[i][i + 1]
        return self.f * sum_s

    def pro_K2(self, hs, he):
        sum_s = 0
        for i in range(self.hs, self.he):
            sum_s += L[i][i + 1]
        return self.m * self.f * sum_s


def read():
    """
    :return: se:可做始终站站号；
    pa_tr:断面客流开行数量
    """
    wb = xlrd.open_workbook('../data/附件1：车站数据.xlsx')
    sheet_st = wb.sheet_by_index(0)
    wb = xlrd.open_workbook('../data/附件2：区间运行时间.xlsx')
    sheet_dur = wb.sheet_by_index(0)
    wb = xlrd.open_workbook('../data/附件3：OD客流数据.xlsx')
    sheet_od = wb.sheet_by_index(0)
    wb = xlrd.open_workbook('../data/附件4：断面客流数据.xlsx')
    sheet_pa = wb.sheet_by_index(0)
    wb = xlrd.open_workbook('../data/附件5：其他数据.xlsx')
    sheet_data = wb.sheet_by_index(0)

    capacity = int(sheet_data.col_values(1)[3])
    # print(capacity)
    launch_str = sheet_data.col_values(1)[0]
    launch_min = launch_str[1:5]
    launch_max = launch_str[5:-1]
    launch = (launch_min, launch_max)
    catch_time = int(sheet_data.col_values(1)[1])
    colpa = sheet_pa.col_values(1)
    pa_tr = np.zeros((32, 32), dtype=np.int_)
    for i in range(1, len(colpa)):
        r_val = colpa[i]
        pa_tr[i][i + 1] = r_val

    bool_station = sheet_st.col_values(1)
    se = []
    for i in range(1, len(bool_station)):
        bs = bool_station[i]
        if bs == "是":
            se.append(i)
    od = np.zeros((31, 31), dtype=np.int_)
    for i in range(1, 31):
        for j in range(i, 31):
            # print(i,j)
            od[i][j] = sheet_od.cell_value(i, j)
    """R , L"""
    oR = np.zeros((31, 31), dtype=np.int_)
    oL = np.zeros((31, 31), dtype=np.float_)
    for i in range(1, 30):
        oR[i][i + 1] = sheet_dur.cell_value(i, 1)
        oL[i][i + 1] = sheet_dur.cell_value(i, 2)

    return se, pa_tr, od, oR, oL


def generate():
    while True:
        f = np.random.randint(1, 30)
        # print(f)
        m = np.random.randint(1, 30)
        i = np.random.randint(1, len(se) - 2)
        j = np.random.randint(i + 1, len(se) - 1)
        if 3 < se[j] - se[i] + 1 <= 24:
            hs = se[i]
            he = se[j]
            if 26 <= (1 + m) * f <= 30:
                gene = Gene(m, f, hs, he)
                gene.constraint()
                if gene.yes:
                    return gene


def initpop(pop_size):
    pop = []
    for i in range(0, pop_size):
        pop.append(generate())
    return pop


def crossover(pop, pc):
    pop_size = len(pop)
    for i in range(0, pop_size - 1):
        if np.random.random() < pc:
            p = pop[i]
            q = pop[i + 1]
            newp = Gene(p.m, p.f, q.hs, q.he)
            newp.constraint()
            if newp.yes:
                pop[i] = newp
            newq = Gene(q.m, q.f, p.hs, p.he)
            newq.constraint()
            if newq.yes:
                pop[i + 1] = newq
    return pop


def mutation(pop, pm):
    pop_size = len(pop)
    for i in range(0, pop_size):
        if np.random.random() < pm:
            gene = generate()
            pop[i] = gene
    return pop


def binary_championship_selection(tmp_pop):
    """
        二元锦标赛排序
        :param tmp_pop:按照快速非支配排序排好的中间种群
        :return:被选中的种群
        """
    size = len(tmp_pop)
    # print("====================================")
    # print("size:",size)
    pop = []
    count = 0
    for i in range(0, size):
        first = np.random.randint(0, size)
        sec = np.random.randint(0, size)
        chosen = first if first <= sec else sec
        # print("chosen==========================", end="")
        # print(tmp_pop[chosen].duration)
        pop.append(tmp_pop[chosen])
        count += 1
        if count == size / 2:
            break
    # print(len(pop))
    # sys.exit(666)
    return pop


def run():
    pop_size = 30
    pc = 0.9
    pm = 0.1
    maxGen = 30
    global od, R, L, se, pa
    se, pa, od, R, L = read()
    # print(pa)
    print(se)
    # print(L)
    # print(od[3][5])
    pop = initpop(pop_size)
    gen = 0
    while gen < maxGen:
        gen = gen + 1
        print("gen: ", gen)
        ancestor = pop.copy()
        for gene in pop:
            gene.fit()
        pop = crossover(pop, pc)
        pop = mutation(pop, pm)
        for gene in pop:
            gene.fit()
        tmp = pop + ancestor
        tmp, paretoRate, paretoNum = NSGA2(tmp)
        pop = binary_championship_selection(tmp)
        """draw"""
        x = [gen] * len(pop)
        y1 = [a.CP for a in pop]
        y2 = [a.CO for a in pop]
        plt.subplot(3, 1, 1)
        plt.xlabel("gen")
        plt.ylabel("CP")
        plt.scatter([gen] * len(y1), y1)
        # 成本散点
        plt.subplot(3, 1, 2)
        plt.xlabel("gen")
        plt.ylabel("CO")
        plt.scatter([gen] * len(y2), y2)
        # 亲子合并中间种群帕累托非支配比例
        plt.subplot(3, 1, 3)
        plt.xlabel("gen")
        plt.ylabel("paretoRate")
        plt.scatter(gen, paretoRate, marker='*')
        plt.grid(True, alpha=0.5, axis='both', linestyle=':')
        """"""
        tmp_res = set(pop)  # 按属性去重，实现ep和hash方法
        global result
        result = result | tmp_res
    plt.savefig("../analysis/res1.png")
    """对最终解集进行排序"""
    result = list(result)
    result, res_paretoRate, paretoNum = NSGA2(result)
    result = result[0:paretoNum]
    for indiv in result:
        print(indiv)
    plt.figure()
    x = [a.CP for a in result]
    y = [a.CO for a in result]
    plt.title("result")
    plt.xlabel("CP")
    plt.ylabel("CO")
    # plt.annotate(txt, (a.resources[0] for a in result))
    plt.scatter(x, y)
    plt.savefig("../analysis/solution1.png")


if __name__ == '__main__':
    run()
    plt.show()
